
##################################################################
##################################################################
## Replication Material
## Simonsen & Widmann: When do Political Parties Moralize?
## A Cross-National Study of the Use of Moral Language in Political Communication on Immigration
## British Journal of Politics
## bakkaer@ps.au.dk
##
## Script 01: Replication of the Main Analysis & Appendix
##################################################################
##################################################################

# Note: The file 000_readme.pdf provides general information to replicate the analysis

# This script was run on the following R version, platform and OS:
# R version 4.3.1 (2023-06-16 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 11 x64 (build 22631)

sessionInfo()


#### Load packages -------------------------------------------------

library(quanteda)
library(tidyverse)
library(tinytable)
library(plm)
library(dotwhisker)
library(broom)
library(summarytools)
library(modelsummary)
library(gridExtra)


### Working Directory ---------------------------------------------

#Set working directory to the downloaded folder
setwd("./files")


### Load Data -----------------------------------------------------

load("./combined_df.Rdata")
load("./ci_combined.Rdata")


### Manuscript Main ------------------------------------------------

##### Table 2 ------------------------------------------------------


# Define file paths for dictionaries and datasets
dictionary_files <- list(
  "G:/My Drive/Work/Moral Dictionary/dictionaries/mrd.yml",
  "G:/My Drive/Work/Moral Dictionary/dictionaries/mrd_da.yml",
  "G:/My Drive/Work/Moral Dictionary/dictionaries/mrd_de.yml",
  "G:/My Drive/Work/Moral Dictionary/dictionaries/mrd_swe.yml",
  "G:/My Drive/Work/Moral Dictionary/dictionaries/mrd_nl.yml"
)

data_files <- list(
  "./cc_english_prepared.Rdata",
  "./cc_danish_prepared.rdata",
  "./cc_german_prepared.Rdata",
  "./cc_swedish_prepared.Rdata",
  "./cc_dutch_prepared.Rdata"
)

languages <- list(
  "english",
  "danish",
  "german",
  "swedish",
  "dutch"
)

stopword_lang <- list(
    "en",
    "da",
    "de",
    "swe",
    "nl"
)

# Variables to analyze
variables <- c("morality")

# Function to apply dictionary and process text
apply_dic <- function(data, dictionary, stopword_lang) {
  corp <- corpus(data)
  toks <- tokens(corp, remove_punct = TRUE, remove_numbers = TRUE, remove_url = TRUE)
  terms_dfm <- dfm(toks)
  toks_removed <- tokens_remove(toks, stopwords(stopword_lang))
  emo_dfm <- dfm(toks_removed)
  dict_dfm_results <- dfm_lookup(emo_dfm, dictionary)
  results_df <- cbind(data, convert(dict_dfm_results, to = "data.frame"))
  results_df$terms_raw <- ntoken(terms_dfm)
  results_df$terms <- ntoken(emo_dfm)
  return(results_df)
}

normalize_and_flag <- function(results) {
  # Explicit normalization
  results$morality.norm <- results[[54]] / results$terms
  
  # Explicit binary flagging
  results$morality_db <- 0
  results$morality_db[results$morality.norm > 0] <- 1
  
  return(results)
}



compute_performance <- function(results, variables) {
  performance <- lapply(variables, function(var) {
    predict <- results[[paste0(var, "_db")]]
    true <- results[[paste0(var, "_binary")]]
    
    # Ensure no division by zero
    retrieved <- sum(predict)
    prec <- if (retrieved > 0) sum(predict & true) / retrieved else 0
    rec <- if (sum(true) > 0) sum(predict & true) / sum(true) else 0
    fscore <- if (prec + rec > 0) 2 * prec * rec / (prec + rec) else 0
    
    data.frame(variable = var, prec = prec, rec = rec, fscore = fscore)
  })
  return(do.call(rbind, performance))
}


# Main loop to process each language
all_performance <- data.frame()

for (lang in 1:5) {
  # Load data
  load(data_files[[lang]])
  
  # Identify the dataset name
  data_name <- ls(pattern = paste0("cc_", languages[[lang]]))
  
  # Assign the dataset to a variable
  data <- get(data_name[1])
  colnames(data)[1] <- "text"  # Rename the first column to "text"
  
  # Load dictionary
  dictionary <- dictionary(file = dictionary_files[[lang]], format = "YAML")
  
  # Apply dictionary and process text
  results <- apply_dic(data, dictionary, tolower(stopword_lang[[lang]]))
  
  # Normalize scores and create binary variables
  results <- normalize_and_flag(results)
  
  # Debugging output
  print(paste("Language:", languages[[lang]]))
  print(head(results$morality.norm))
  print(table(results$morality_db))
  
  # Compute performance
  df_performance <- compute_performance(results, variables)
  df_performance <- df_performance[1, ]  # Keep only one row
  df_performance$language <- languages[[lang]]
  
  # Combine performance results
  all_performance <- rbind(all_performance, df_performance)
}

# Final output
print(all_performance)

# Create a flextable for better formatting and save as Word document
all_performance %>%
  tt(rownames = TRUE) %>%   # Convert to formatted table with row names
  save_tt("./output/table2.docx", overwrite = TRUE)



##### Figure 1 -----------------------------------------------------

# Fit a panel data model with fixed effects
model1 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)

# Summarize the model results
summary(model1)

# Tidy the model output and relabel predictors for better readability
tidy_model1 <- tidy(model1) %>%
  relabel_predictors(c(
    immipos_year_scaled = "Immigration Position",
    extremism_year_scaled2 = "Extremity",
    gov1 = "Government",
    polarization_year_scaled = "Polarization",
    salience_year_scaled = "Salience",
    election = "Election Year"
  ))

# Calculate 95% confidence intervals for model coefficients
fit_cis_95 <- confint(model1, level = 0.95) %>% 
  as.data.frame() %>%
  rename(conf.low_95 = `2.5 %`, conf.high_95 = `97.5 %`)

# Calculate 90% confidence intervals for model coefficients
fit_cis_90 <- confint(model1, level = 0.90) %>% 
  as.data.frame() %>%
  rename(conf.low_90 = `5 %`, conf.high_90 = `95 %`)

# Combine tidy results with confidence intervals and clean column names
tidy_model1 <- tidy_model1 %>%
  bind_cols(fit_cis_95, fit_cis_90) %>%
  rename(
    Variable = term,
    Coefficient = estimate,
    SE = std.error
  ) %>%
  filter(Variable != "(Intercept)")  # Remove the intercept row for clarity

# Create a plot of the coefficients with confidence intervals
tidy_model1 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>%  # Order variables for display
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(
    yintercept = 0, 
    colour = gray(0.5), 
    linetype = "dashed"
  ) +  # Add a reference line at 0
  geom_point(
    aes(x = Variable, y = Coefficient), 
    color = "coral"
  ) +  # Add points for coefficients
  geom_linerange(
    aes(x = Variable, ymin = conf.low_95, ymax = conf.high_95),
    color = "coral", 
    linewidth = 0.5
  ) +  # Add 95% confidence intervals
  xlab("") + 
  ylab("") + 
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +  # Use a clean theme
  theme(
    text = element_text(size = 20),
    legend.position = "none",
    plot.title = element_text(size = 15, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 1),
    plot.caption = element_text(size = 10, hjust = 0)
  ) +
  coord_flip()  # Flip coordinates for better readability


# saving Figure 1
ggsave(
  filename = "./output/figure1.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


### Supplementary Material ----------------------------------------


#### Appendix D ---------------------------------------------------

##### Table D1 ---------------------
# Generate descriptive statistics for selected variables
table1 <- descr(
  combined_df[c("morality_year", "extremism2", "immipos", "polarization", "salience_year")], 
  stats = c("mean", "sd", "min", "med", "max")
) %>% 
  as.data.frame()

# Rename rows for better readability
colnames(table1) <- c(
  "Extremity",
  "Immigration Position",
  "Moral Rhetoric",
  "Polarization",
  "Salience"
)

# Format all numerical values to 2 decimals
table1 <- table1 %>%
  mutate(across(everything(), ~ round(.x, 2)))  # Applies rounding to all columns

# Create a flextable for better formatting and save as Word document
table1 %>%
  tt(rownames = TRUE) %>%   # Convert to formatted table with row names
  save_tt("./output/tableD1.docx", overwrite = TRUE)


##### Figure D1 ---------------------

ggplot(combined_df, aes(x=morality_year)) + 
  geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                 colour="black", fill="white") +
  geom_density(alpha=.5, fill="coral") +
  xlab("Moral Rhetoric") +
  theme_bw()

# saving figure 
ggsave(
  filename = "./output/figureD1.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)

##### Figure D2 ---------------------

ggplot(combined_df, aes(x=immipos)) + 
  geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                 colour="black", fill="white") +
  geom_density(alpha=.5, fill="coral") +
  xlab("Immigration Position") +
  theme_bw()

# saving figure 
ggsave(
  filename = "./output/figureD2.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure D3 ---------------------

ggplot(combined_df, aes(x=extremism2)) + 
  geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                 colour="black", fill="white") +
  geom_density(alpha=.5, fill="coral") +
  xlab("Extremity") +
  theme_bw()

# saving figure 
ggsave(
  filename = "./output/figureD3.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure D4 ---------------------

ggplot(combined_df, aes(x=polarization)) + 
  geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                 colour="black", fill="white") +
  geom_density(alpha=.5, fill="coral") +
  xlab("Polarization") +
  theme_bw()

# saving figure 
ggsave(
  filename = "./output/figureD4.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)



##### Figure D5 ---------------------

ggplot(combined_df, aes(x=salience_year)) + 
  geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                 colour="black", fill="white") +
  geom_density(alpha=.5, fill="coral") +
  xlab("Salience") +
  theme_bw()


# saving figure 
ggsave(
  filename = "./output/figureD5.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)



##### Figure D6 ---------------------

# Create a named vector for custom labels
country_labels <- c(
  "austria" = "Austria",
  "canada" = "Canada",
  "denmark" = "Denmark",
  "germany" = "Germany",
  "netherlands" = "Netherlands",
  "sweden" = "Sweden",
  "uk" = "UK",
  "us" = "US")


ggplot(combined_df, aes(x = morality_year)) +
  geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
  geom_density(alpha = 0.5, fill = "coral") +
  facet_wrap(~ country, labeller = labeller(country = country_labels)) +
  xlab("Moral Rhetoric")

# saving figure 
ggsave(
  filename = "./output/figureD6.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure D7 ---------------------


combined_df <- combined_df %>%
  mutate(year = as.numeric(as.character(year)),  # Convert to numeric if it's not
         decade = floor(year / 10) * 10)  # Now create the 'decade' variable


# Now, plot the histogram using the new 'decade' variable
ggplot(combined_df, aes(x=morality_year)) +
  geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                 colour="black", fill="white") +
  geom_density(alpha=.5, fill="coral") +
  facet_wrap(~decade) +   
  xlab("Decade") +
  theme_bw() 

# saving figure 
ggsave(
  filename = "./output/figureD7.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure D8 -------------------------

combined_subset <- combined_df[!is.na(combined_df$partytype),]

ggplot(combined_subset, aes(x = morality_year)) +
  geom_histogram(aes(y = ..density..), colour = "black", fill = "white") +
  geom_density(alpha = 0.5, fill = "coral") +
  facet_wrap(~ partytype) +
  xlab("Moral Rhetoric")

# saving figure 
ggsave(
  filename = "./output/figureD8.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure D9 --------------------------

# Aggregate results by category
ci_combined2 <- ci_combined %>%
  group_by(category) %>%
  summarise(mean_morality = mean(mean_morality, na.rm = TRUE),
            lower_ci = mean(lower_ci, na.rm = TRUE),
            upper_ci = mean(upper_ci, na.rm = TRUE))

# Plot the overall means with error bars using ggplot2
ggplot(ci_combined2, aes(x = category, y = mean_morality, fill = category)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7, alpha = 0.7) +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.25) +
  labs(title = "Comparison of Moral Rhetoric",
       x = "Speech Category",
       y = "Average Moral Rhetoric") +
  scale_fill_manual(values = c("Immigration" = "#0072B2", "Non-Immigration" = "coral"),
                    labels = c("Immigration", "Non-Immigration")) +
  theme_minimal()

# saving figure 
ggsave(
  filename = "./output/figureD9.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


#### Appendix E ---------------------------------------------------


##### Table E1 ---------------------


# Fit a panel data model with fixed effects
model1 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)


# Fit a panel data model with fixed effects
model2 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + weighted_polarization_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)


# Fit a panel data model with fixed effects
model3 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + pol_cult_vdem + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)

# Fit a panel data model with fixed effects
model4 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "party",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)

combined_df <- combined_df |> 
  group_by(party) |> 
  arrange(year) |> 
  mutate(lag_morality = dplyr::lag(morality_year_scaled)) |> 
  ungroup()

# Fit a panel data model with fixed effects
model5 <- plm(
  lag_morality ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)

# with splitting sample
combined_df_tvearly <- combined_df %>%
  filter(televised == 0)

combined_df_tvlate <- combined_df %>%
  filter(televised == 1)


# Fit a panel data model with fixed effects
model6 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df_tvearly  # Data frame used for analysis
)

# Fit a panel data model with fixed effects
model7 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df_tvlate  # Data frame used for analysis
)



combined_df_smearly <- combined_df %>%
  filter(socialmedia == 0)

combined_df_smlate <- combined_df %>%
  filter(socialmedia == 1)


# Fit a panel data model with fixed effects
model8 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df_smearly  # Data frame used for analysis
)

# Fit a panel data model with fixed effects
model9 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df_smlate  # Data frame used for analysis
)



# Fit a panel data model with fixed effects
model10 <- plm(
  morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 +
    gov1 + polarization_year_scaled + salience_year_scaled + election + sentiment_year_scaled,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)

model11 <- plm(
  morality_year_scaled ~ gov1 + polarization_year_scaled + salience_year_scaled + election + partytype,
  index = "country",  # Specifies the panel structure
  model = "within",   # Fixed effects model
  data = combined_df  # Data frame used for analysis
)



models <- list(model1, model2, model3, model4, model5, model6, model7, model8, model9, model10, model11)
modelsummary(models, stars = TRUE)
modelsummary(models, stars = TRUE, output = "./output/tableE1.docx")


##### Figure E1 ---------------------


tidy2 <- tidy(model2) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       weighted_polarization_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model2, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model2, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy2 <- bind_cols(tidy2, 
                     fit_cis_95, 
                     fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


tidy2 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()

# saving figure 
ggsave(
  filename = "./output/figureE1.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure E2 ---------------------

tidy3 <- tidy(model3) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       pol_cult_vdem = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model3, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model3, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy3 <- bind_cols(tidy3, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


tidy3 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()

# saving figure 
ggsave(
  filename = "./output/figureE2.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure E3 ---------------------


tidy4 <- tidy(model4) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model4, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model4, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy4 <- bind_cols(tidy4, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


tidy4 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()

# saving figure 
ggsave(
  filename = "./output/figureE3.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)

##### Figure E4 ---------------------


tidy5 <- tidy(model5) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model5, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model5, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy5 <- bind_cols(tidy5, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


tidy5 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()

# saving figure 
ggsave(
  filename = "./output/figureE4.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)

##### Figure E5 ---------------------


tidy6 <- tidy(model6) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model6, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model6, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy6 <- bind_cols(tidy6, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


p1 <- tidy6 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()





tidy7 <- tidy(model7) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model7, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model7, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy7 <- bind_cols(tidy7, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


p2 <- tidy7 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()



arranged_plot <- grid.arrange(p1, p2, ncol = 2)  # Arrange plots side by side


# saving figure 
ggsave(
  plot = arranged_plot,
  filename = "./output/figureE5.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)



##### Figure E6 ---------------------


tidy8 <- tidy(model8) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model8, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model8, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy8 <- bind_cols(tidy8, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


p1 <- tidy8 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()





tidy9 <- tidy(model9) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year"))

fit_cis_95 <- confint(model9, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model9, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy9 <- bind_cols(tidy9, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


p2 <- tidy9 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()



arranged_plot <- grid.arrange(p1, p2, ncol = 2)  # Arrange plots side by side


# saving figure 
ggsave(
  plot = arranged_plot,
  filename = "./output/figureE6.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure E7 ---------------------


tidy10 <- tidy(model10) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year",
                       sentiment_year_scaled = "Sentiment"))

fit_cis_95 <- confint(model10, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model10, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy10 <- bind_cols(tidy10, 
                   fit_cis_95, 
                   fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


tidy10 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()

# saving figure 
ggsave(
  filename = "./output/figureE7.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure E8 ---------------------


tidy11 <- tidy(model11) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year",
                       sentiment_year_scaled = "Sentiment",
                       partytypemoderate = "Moderates",
                       partytypepro = "Pro"))

fit_cis_95 <- confint(model11, level = 0.95) %>% 
  data.frame() %>%
  rename("conf.low_95" = "X2.5..",
         "conf.high_95" = "X97.5..")
fit_cis_90 <- confint(model11, level = 0.90) %>% 
  data.frame() %>%
  rename("conf.low_90" = "X5..",
         "conf.high_90" = "X95..")

tidy11 <- bind_cols(tidy11, 
                    fit_cis_95, 
                    fit_cis_90) %>%
  rename(Variable = term,
         Coefficient = estimate,
         SE = std.error) %>%
  filter(Variable != "(Intercept)")


tidy11 %>%
  mutate(Variable = fct_reorder(Variable, desc(Variable))) %>% 
  ggplot(aes(x = Variable, y = Coefficient)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = Variable, 
                 y = Coefficient), color = "coral") + 
  geom_linerange(aes(x = Variable, 
                     ymin = conf.low_95,
                     ymax = conf.high_95, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0)) +
  coord_flip()

# saving figure 
ggsave(
  filename = "./output/figureE8.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


#### Appendix F ---------------------------------------------------

##### Figure F1 ----------------------

fiti1 <- plm(morality_year_scaled ~ extremism_year_scaled2  + immipos_year_scaled*polarization_year_scaled + salience_year_scaled + election + gov1, index = "country", model = "within", data = combined_df)

model_tidy <- tidy(fiti1) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year",
                       'immipos_year_scaled:polarization_year_scaled' = "Immigration Position*Polarization"))

model_tidy %>%
  mutate(term = fct_reorder(term, desc(term))) %>% 
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = estimate, 
                 y = term), color = "coral") + 
  geom_linerange(aes(x = estimate, 
                     xmin = estimate - std.error * 1.96,
                     xmax = estimate + std.error * 1.96, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0))

# saving figure 
ggsave(
  filename = "./output/figureF1.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)



##### Figure F2 ----------------------


fiti2 <- plm(morality_year_scaled ~ extremism_year_scaled2*polarization_year_scaled  + immipos_year_scaled + salience_year_scaled + election + gov1, index = "country", model = "within", data = combined_df)

model_tidy <- tidy(fiti2) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year",
                       'extremism_year_scaled2:polarization_year_scaled' = "Extremity*Polarization"))

model_tidy %>%
  mutate(term = fct_reorder(term, desc(term))) %>% 
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = estimate, 
                 y = term), color = "coral") + 
  geom_linerange(aes(x = estimate, 
                     xmin = estimate - std.error * 1.96,
                     xmax = estimate + std.error * 1.96, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0))

# saving figure 
ggsave(
  filename = "./output/figureF2.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Figure F3 ----------------------


fiti3 <- plm(morality_year_scaled ~ extremism_year_scaled2  + immipos_year_scaled + salience_year_scaled + election + gov1*polarization_year_scaled, index = "country", model = "within", data = combined_df)

model_tidy <- tidy(fiti3) %>%
  relabel_predictors(c(immipos_year_scaled = "Immigration Position",
                       extremism_year_scaled2 = "Extremity",
                       gov1 = "Government",
                       polarization_year_scaled = "Polarization",
                       salience_year_scaled = "Salience",
                       election = "Election Year",
                       'gov1:polarization_year_scaled' = "Government*Polarization"))

model_tidy %>%
  mutate(term = fct_reorder(term, desc(term))) %>% 
  ggplot(aes(x = estimate, y = term)) +
  geom_vline(xintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_point(aes(x = estimate, 
                 y = term), color = "coral") + 
  geom_linerange(aes(x = estimate, 
                     xmin = estimate - std.error * 1.96,
                     xmax = estimate + std.error * 1.96, color = "coral"),
                 lwd = 1/2) + 
  xlab("") + ylab("") +
  ggtitle("Predicting Moral Rhetoric") +
  theme_bw() +
  theme(text = element_text(size = 20), legend.position="none", plot.title = element_text(size = 15, hjust = 0.5), axis.title.x=element_text(size = 15, vjust=1), 
        plot.caption = element_text(size = 10, hjust = 0))


# saving figure 
ggsave(
  filename = "./output/figureF3.png",   # Specify the file name and format
  device = "png",            # File format
  dpi = 600,                 # Resolution in dots per inch
  width = 10,                 # Width of the plot in inches (adjust as needed)
  height = 10,                # Height of the plot in inches (adjust as needed)
  units = "in"               # Units for width and height
)


##### Table F1 ------------------

fitF1 <- plm(morality_year_scaled ~ immipos_year_scaled + extremism_year_scaled2 + gov1 + election, index = "country", model = "within", data = combined_df)
modelsummary(fitF1, stars = TRUE)
modelsummary(fitF1, stars = TRUE, output = "./output/tableF1.docx")


#####################################################################################################################
############ END OF SCRIPT ##########################################################################################
#####################################################################################################################
